home *** CD-ROM | disk | FTP | other *** search
/ Atari Mega Archive 1 / Atari Mega Archive - Volume 1.iso / mint / lib / mntlib44.zoo / mntlib / _muldf3.cpp < prev    next >
Text File  |  1993-06-17  |  6KB  |  265 lines
  1. | mjr: not needed on the TT
  2.  
  3. #ifndef    __M68881__
  4.     .text
  5.     .even
  6.     .globl    __muldf3, ___muldf3
  7.  
  8. __muldf3:
  9. ___muldf3:
  10.  
  11. # ifdef    sfp004
  12.  
  13. | double precision floating point stuff for Atari-gcc using the SFP004
  14. | developed with gas
  15. |
  16. | double precision multiplication
  17. |
  18. | M. Ritzert (mjr at dfg.dbp.de)
  19. |
  20. | 4.10.1990
  21. |
  22. | no NAN checking implemented since the 68881 treats this situation "correctly",
  23. | i.e. according to IEEE
  24.  
  25. | addresses of the 68881 data port. This choice is fastest when much data is
  26. | transferred between the two processors.
  27.  
  28. comm =     -6
  29. resp =    -16
  30. zahl =      0
  31.  
  32. | waiting loop ...
  33. |
  34. | wait:
  35. | ww:    cmpiw    #0x8900,a0@(resp)
  36. |     beq    ww
  37. | is coded directly by
  38. |    .long    0x0c688900, 0xfff067f8
  39.  
  40.     lea    0xfffffa50:w,a0
  41.     movew    #0x5400,a0@(comm)    | load first argument to fp0
  42.     cmpiw    #0x8900,a0@(resp)    | check
  43.     movel    a7@(4),a0@
  44.     movel    a7@(8),a0@
  45.     movew    #0x5423,a0@(comm)
  46.     .long    0x0c688900, 0xfff067f8
  47.     movel    a7@(12),a0@
  48.     movel    a7@(16),a0@
  49.     movew    #0x7400,a0@(comm)    | result to d0/d1
  50.     .long    0x0c688900, 0xfff067f8
  51.     movel    a0@,d0
  52.     movel    a0@,d1
  53.     rts
  54.  
  55. # else    sfp004
  56.  
  57. | double floating point multiplication routine
  58. |
  59. | written by Kai-Uwe Bloem (I5110401@dbstu1.bitnet).
  60. | Based on a 80x86 floating point packet from comp.os.minix, written by P.Housel
  61. | Revision 1.2.4 michal 05-93 (ntomczak@vm.ucs.ualberta.ca)
  62. |   + ensure that Inf * NaN == NaN * Inf == NaN
  63. |     and 0 * Inf = Inf * 0 = NaN
  64. |
  65. | Revision 1.2.3 michal 05-93 (ntomczak@vm.ucs.ualberta.ca)
  66. |   + code smoothing
  67. |
  68. | patched by Olaf Flebbe (flebbe@tat.physik.uni-tuebingen.de)
  69. |
  70. | Revision 1.2.2 olaf 05-93:
  71. |   + fixed a bug for signed bug for 0.
  72. |
  73. | Revision 1.2.1 olaf 12-92:
  74. |   + added support for NaN and Infinites
  75. |   + added support for -0
  76. |
  77. | Revision 1.2, kub 01-90 :
  78. | added support for denormalized numbers
  79. |
  80. | Revision 1.1, kub 12-89 :
  81. | Ported over to 68k assembler
  82. |
  83. | Revision 1.0:
  84. | original 8088 code from P.S.Housel
  85.  
  86. BIAS8    =    0x3FF-1
  87.  
  88.     lea    sp@(4),a0
  89.     moveml    d2-d7,sp@-
  90.     moveml    a0@,d4-d5/d6-d7 | d4-d5 = v, d6-d7 = u
  91.  
  92.     movel    #0x0fffff,d3
  93.     movel    d6,d0        | d0 = u.exp
  94.     andl    d3,d6        | remove exponent from u.mantissa
  95.     swap    d0
  96.     movew    d0,d2        | d2 = u.sign
  97.  
  98.     movel    d4,d1        | d1 = v.exp
  99.     andl    d3,d4        | remove exponent from v.mantissa
  100.     swap    d1
  101.     eorw    d1,d2        | d2 = u.sign ^ v.sign (in bit 15)
  102.     
  103.     moveq    #15,d3
  104.     bclr    d3,d0        | kill sign bit
  105.     bclr    d3,d1        | kill sign bit
  106.     tstl    d0        | test if one of factors is 0
  107.     bne    0f        | if not the first one then maybe the second
  108.     tstl    d7
  109.     beq    1f
  110. 0:    tstl    d1
  111.     bne    1f
  112.     tstl    d5
  113. 1:    seq    d2        | 'one of factors is 0' flag in the lowest byte
  114.     lsrw    #4,d0        | keep in d0, d1 exponents only
  115.     lsrw    #4,d1
  116. |
  117. | Testing for NaN and Infinities
  118. |
  119.     movew    #0x7ff,d3
  120.     cmpw    d3,d0
  121.     beq    0f
  122.     cmpw    d3,d1
  123.     bne    nospec
  124.     bra    1f
  125. |    first operand is special
  126. |    Nan?
  127. 0:    orl    d7,d6
  128.     bne    retnan
  129. 1:    tstb    d2        | 0 times special or special times 0 ?
  130.     bne    retnan        | yes -> NaN
  131.     cmpw    d3,d1        | is the other special?
  132.     beq    2f        | maybe it is NaN
  133. |
  134. |    Return Infinity with correct sign
  135. |
  136. retinf:    moveq    #0,d1
  137.     movel    #0xffe00000,d0    | we will return #0xfff00000 or #0x7ff00000
  138.     lslw    #1,d2
  139.     roxrl   #1,d0        | shift in high bit as given by d2
  140. return:    moveml    sp@+,d2-d7
  141.     rts
  142. |
  143. | v is special
  144. |
  145. 2:    orl    d5,d4        | is v NaN?
  146.     beq    retinf        | if not then we have (not-NaN * Inf)
  147. |
  148. |    Return NaN
  149. |
  150. retnan: moveql    #-1,d1
  151.     movel    d1,d0
  152.     lsrl    #1,d0        | 0x7fffffff -> d0
  153.     bra    return
  154. |
  155. | end of NaN and Inf.
  156. |
  157. nospec:    tstb    d2        | not needed - but we can waste two instr.
  158.     bne    retzz        | return signed 0 if one of factors is 0
  159.     lea    sp@(-16),sp    | multiplication accumulator
  160.  
  161.     moveq    #20,d3
  162.     bset    d3,d6        | restore implied leading "1"
  163.     tstw    d0        | check for zero exponent - no leading "1"
  164.     bne    1f
  165.     bclr    d3,d6        | remove it
  166.     addqw    #1,d0        | "normalize" exponent
  167. 1:    movel    d6,d3
  168.     orl    d7,d3
  169.     beq    retz        | multiplying by zero
  170.  
  171.     moveq    #20,d3
  172.     bset    d3,d4        | restore implied leading "1"
  173.     tstw    d1        | check for zero exponent - no leading "1"
  174.     bne    1f
  175.     bclr    d3,d4        | remove it
  176.     addqw    #1,d1        | "normalize" exponent
  177. 1:    movel    d4,d3
  178.     orl    d5,d3
  179.     beq    retz        | multiplying by zero
  180.  
  181.     addw    d1,d0        | add exponents,
  182.     subw    #BIAS8+16-11,d0    | remove excess bias, acnt for repositioning
  183.  
  184.     lea    sp@,a1        | initialize 128-bit product to zero
  185.     moveq    #0,d3
  186.     movel    d3,a1@+
  187.     movel    d3,a1@+
  188.     movel    d3,a1@+
  189.     movel    d3,a1@
  190.     subqw    #4,a1        | an address of sp@(8) in a1
  191.  
  192. | see Knuth, Seminumerical Algorithms, section 4.3. algorithm M
  193.  
  194.     swap    d2
  195.     movew    #4-1,d2
  196. 1:
  197.     movew    d5,d3
  198.     mulu    d7,d3        | mulitply with bigit from multiplier
  199.     addl    d3,a1@(4)    | store into result
  200.     movew    d4,d3
  201.     mulu    d7,d3
  202.     movel    a1@,d1        | add to result
  203.     addxl    d3,d1
  204.     movel    d1,a1@
  205.     roxlw    a1@-        | rotate carry in
  206.  
  207.     movel    d5,d3
  208.     swap    d3
  209.     mulu    d7,d3
  210.     addl    d3,a1@(4)    | add to result
  211.     movel    d4,d3
  212.     swap    d3
  213.     mulu    d7,d3
  214.     movel    a1@,d1        | add to result
  215.     addxl    d3,d1
  216.     movel    d1,a1@
  217.  
  218.     movew    d6,d7
  219.     swap    d6
  220.     swap    d7
  221.     dbra    d2,1b
  222.  
  223.     swap    d2        | [TOP 16 BITS SHOULD BE ZERO !]
  224.  
  225.     moveml    sp@(2),d4-d7    | get the 112 valid bits
  226.     clrw    d7        | (pad to 128)
  227.     movel    #0x0000ffff,d3
  228. 2:
  229.     cmpl    d3,d4        | multiply (shift) until
  230.     bhi    3f        |  1 in upper 16 result bits
  231.     cmpw    #9,d0        | give up for denormalized numbers
  232.     ble    3f
  233.     swap    d4        | (we''re getting here only when multiplying
  234.     swap    d5        |  with a denormalized number; there''s an
  235.     swap    d6        |  eventual loss of 4 bits in the rounding
  236.     swap    d7        |  byte -- what a pity 8-)
  237.     movew    d5,d4
  238.     movew    d6,d5
  239.     movew    d7,d6
  240.     clrw    d7
  241.     subqw    #8,d0        | decrement exponent
  242.     subqw    #8,d0
  243.     bra    2b
  244. 3:
  245.     movel    d6,d1        | get rounding bits
  246.     roll    #8,d1
  247.     movel    d1,d3        | see if sticky bit should be set
  248.     orl    d7,d3        | (lower 16 bits of d7 are guaranteed to be 0)
  249.     andl    #0xffffff00,d3
  250.     beq    4f
  251.     orb    #1,d1        | set "sticky bit" if any low-order set
  252. 4:    lea    sp@(16),sp    | remove accumulator from stack
  253.     jmp    norm_df        | (result in d4/d5)
  254. |                | norm_df does not return here
  255. retz:    lea    sp@(16),sp    | drop accumulator space
  256. retzz:    moveq    #0,d0        | save zero as result
  257.     movel    d0,d1
  258.     lslw    #1,d2        | fill X bit
  259.     roxrl    #1,d0        | set high bit of d0 accordingly
  260.     moveml    sp@+,d2-d7
  261.     rts            | no normalizing neccessary
  262.  
  263. # endif    sfp004
  264. #endif    __M68881__
  265.